#include <iostream>
#include<sys/time.h>
#include<pthread.h>
#include<semaphore.h>
#include<xmmintrin.h>
#include<immintrin.h>
#include<x86intrin.h>
using namespace std;
#define N 4800
#define num_thread 8
int n=100;
float m[N][N],init_m[N][N];
void init(int n){
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)
            init_m[i][j]=0;
        init_m[i][i]=1.0;
        for(int j=i+1;j<n;j++)
            init_m[i][j]=rand()%1000+1;
    }
    for(int k=0;k<n;k++)
        for(int i=k+1;i<n;i++)
            for(int j=0;j<n;j++)
                init_m[i][j]=int((init_m[i][j]+init_m[k][j]))%1000+1.0;
}
void reset(int n){
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            m[i][j]=init_m[i][j];
}
void serial(int n){ //串行
    struct timeval h, t;
    float time = 0.0;
    gettimeofday(&h, NULL);
      for(int k=0;k<n;k++){
          for(int j=k+1;j<n;j++)
              m[k][j]/=m[k][k];
          m[k][k]=1.0;
          for(int i=k+1;i<n;i++){
              for(int j=k+1;j<n;j++)
                  m[i][j]-=m[i][k]*m[k][j];
              m[i][k]=0;
          }
      }
    gettimeofday(&t, NULL);
    time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"serial: "<<time<<" ms"<<endl;
    return;
}


typedef struct {
    int k;//消去的轮次
    int t_id;//线程id
}threadParam_t;
sem_t sem_main; //信号量定义
sem_t sem_workerstart[num_thread];
sem_t sem_workerend[num_thread];

//静态线程行分配SSE
void* threadFunc1(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;
    for(int k = 0; k < n; k++) {
        sem_wait(&sem_workerstart[t_id]);//阻塞，等待主线程完成除法操作
        //循环划分任务
        for(int i = k + t_id + 1; i < n; i += num_thread) {
            float* aik = new float[4];
            *aik = *(aik+1) = *(aik+2) = *(aik+3) = m[i][k];
            __m128 vaik = _mm_load_ps(aik);//将四个单精度浮点数从内存加载到向量寄存器
            delete[] aik;
            for(int j = k + 1; j + 3 < n; j += 4) {
                __m128 vaij = _mm_loadu_ps(&m[i][j]);
                __m128 vakj = _mm_loadu_ps(&m[k][j]);
                __m128 vx = _mm_mul_ps(vaik, vakj);
                vaij = _mm_sub_ps(vaij, vx);//A[i][j] = A[i][j] - A[i][k]*A[k][j];
                _mm_storeu_ps(&m[i][j], vaij);//存储到内存
            }
            for(int j = n - n%4; j < n; j++) {
                m[i][j] = m[i][j] - m[i][k]*m[k][j];
            }//结尾几个元素串行计算
            m[i][k] = 0;
        }
        sem_post(&sem_main);//唤醒主线程
        sem_wait(&sem_workerend[t_id]);//阻塞，等待主线程唤醒进入下一轮
    }
    pthread_exit(nullptr);
}
//静态线程行分配AVX
void* threadFunc2(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;
    for(int k = 0; k < n; k++) {
        sem_wait(&sem_workerstart[t_id]);//阻塞，等待主线程完成除法操作
        //循环划分任务
        for(int i = k + t_id + 1; i < n; i += num_thread) {
            float* aik = new float[8];
            *aik = *(aik+1) = *(aik+2) = *(aik+3) = *(aik+4) = *(aik+5) = *(aik+6) = *(aik+7) = m[i][k];
            __m256 vaik = __builtin_ia32_loadups256(aik);//将八个单精度浮点数从内存加载到向量寄存器
            delete[] aik;
            for(int j = k + 1; j + 7 < n; j += 8) {
                __m256 vaij = __builtin_ia32_loadups256(&m[i][j]);
                __m256 vakj = __builtin_ia32_loadups256(&m[k][j]);
                __m256 vx = __builtin_ia32_mulps256(vaik, vakj);
                vaij = __builtin_ia32_subps256(vaij, vx);//A[i][j] = A[i][j] - A[i][k]*A[k][j];
                __builtin_ia32_storeups256(&m[i][j], vaij);//存储到内存
            }
            for(int j = n - n%8; j < n; j++) {
                m[i][j] = m[i][j] - m[i][k]*m[k][j];
            }//结尾几个元素串行计算
            m[i][k] = 0;
        }
        sem_post(&sem_main);//唤醒主线程
        sem_wait(&sem_workerend[t_id]);//阻塞，等待主线程唤醒进入下一轮
    }
    pthread_exit(nullptr);
}

void pth1() {
    struct timeval h, t;
    float time = 0.0;
    cout<<"Pthread3:(N = "<<n<<" and P = "<<num_thread<<")\n";
        init(n);
        gettimeofday(&h, NULL);
        //初始化信号量
        sem_init(&sem_main, 0, 0);
        for(int i = 0; i < num_thread; i++) {
            sem_init(&sem_workerstart[i], 0, 0);
            sem_init(&sem_workerend[i], 0, 0);
        }
        //创建线程
        pthread_t handles[num_thread];
        threadParam_t param[num_thread];
        for(int t_id = 0; t_id < num_thread; t_id++) {
            param[t_id].t_id = t_id;
            pthread_create(handles + t_id, nullptr, threadFunc1, param + t_id);
        }
        for(int k = 0; k < n; k++) {
            //主线程做除法操作
            float* akk = new float[4];
            *akk = *(akk+1) = *(akk+2) = *(akk+3) = m[k][k];
            __m128 vt = _mm_load_ps(akk);//将四个单精度浮点数从内存加载到向量寄存器
            delete[] akk;
            for(int j = k + 1; j + 3 < n; j += 4) {
                __m128 va = _mm_loadu_ps(&m[k][j]);
                va = _mm_div_ps(va, vt);//A[k][j] = A[k][j]/A[k][k];
                _mm_storeu_ps(&m[k][j], va);//储存到内存
            }
            for(int j = n - n%4; j < n; j++) {
                m[k][j] = m[k][j]/m[k][k];
            }//结尾几个元素串行计算
            m[k][k] = 1.0;
            //开始唤醒工作线程
            for(int t_id = 0; t_id < num_thread; t_id++) {
                sem_post(&sem_workerstart[t_id]);
            }
            //主线程睡眠
            for(int t_id = 0; t_id < num_thread; t_id++) {
                sem_wait(&sem_main);
            }
            //主线程再次唤醒工作线程
            for(int t_id = 0; t_id < num_thread; t_id++) {
                sem_post(&sem_workerend[t_id]);
            }
        }
        for(int t_id = 0; t_id < num_thread; t_id++) {
            pthread_join(handles[t_id], nullptr);
        }
        sem_destroy(&sem_main);
        for(int i = 0; i < num_thread; i++) {
            sem_destroy(&sem_workerstart[i]);
            sem_destroy(&sem_workerend[i]);
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"TIME: "<<time<<" ms"<<endl;
}
void pth2() {
    struct timeval h, t;
    float time = 0.0;
    cout<<"Pthread4:(N = "<<n<<" and P = "<<num_thread<<")\n";
        init(n);
        gettimeofday(&h, NULL);
        //初始化信号量
        sem_init(&sem_main, 0, 0);
        for(int i = 0; i < num_thread; i++) {
            sem_init(&sem_workerstart[i], 0, 0);
            sem_init(&sem_workerend[i], 0, 0);
        }
        //创建线程
        pthread_t handles[num_thread];
        threadParam_t param[num_thread];
        for(int t_id = 0; t_id < num_thread; t_id++) {
            param[t_id].t_id = t_id;
            pthread_create(handles + t_id, nullptr, threadFunc2, param + t_id);
        }
        for(int k = 0; k < n; k++) {
            //主线程做除法操作
            float* akk = new float[8];
            *akk = *(akk+1) = *(akk+2) = *(akk+3) = *(akk+4) = *(akk+5) = *(akk+6) = *(akk+7) = m[k][k];
            __m256 vt = __builtin_ia32_loadups256(akk);//将八个单精度浮点数从内存加载到向量寄存器
            delete[] akk;
            for(int j = k + 1; j + 7 < n; j += 8) {
                __m256 va = __builtin_ia32_loadups256(&m[k][j]);
                va = _mm256_div_ps(va, vt);//A[k][j] = A[k][j]/A[k][k];
                __builtin_ia32_storeups256(&m[k][j], va);//储存到内存
            }
            for(int j = n - n%8; j < n; j++) {
                m[k][j] = m[k][j]/m[k][k];
            }//结尾几个元素串行计算
            m[k][k] = 1.0;
            //开始唤醒工作线程
            for(int t_id = 0; t_id < num_thread; t_id++) {
                sem_post(&sem_workerstart[t_id]);
            }
            //主线程睡眠
            for(int t_id = 0; t_id < num_thread; t_id++) {
                sem_wait(&sem_main);
            }
            //主线程再次唤醒工作线程
            for(int t_id = 0; t_id < num_thread; t_id++) {
                sem_post(&sem_workerend[t_id]);
            }
        }
        for(int t_id = 0; t_id < num_thread; t_id++) {
            pthread_join(handles[t_id], nullptr);
        }
        sem_destroy(&sem_main);
        for(int i = 0; i < num_thread; i++) {
            sem_destroy(&sem_workerstart[i]);
            sem_destroy(&sem_workerend[i]);
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"TIME: "<<time<<" ms"<<endl;
}
int main(){
    cout<<"n="<<n<<endl;
    init(n);
    reset(n);
    serial(n);
    reset(n);
    pth1();
    reset(n);
    pth2();
    n=200;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pth1();
    reset(n);
    pth2();
    n=400;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pth1();
    reset(n);
    pth2();
    n=800;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pth1();
    reset(n);
    pth2();
    n=1600;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pth1();
    reset(n);
    pth2();
    n=3200;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pth1();
    reset(n);
    pth2();
    n=4800;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pth1();
    reset(n);
    pth2();
    return 0;
}